Basic map
Add elements to map
plot(wrld_simpl,xlim=c(115,128) ,ylim=c(19.5,27.5),col='#D2B48C',bg='lightblue') # TW map
coords <- matrix(c(121.537290,120.265541, 25.021335, 22.626524),ncol=2) # NTU and SYS univ.
coords <- coordinates(coords) # assign values as spatial coordinates
spoints <- SpatialPoints(coords) # create SpatialPoints
df <- data.frame(location=c("NTU","SYS")) # create a dataframe
spointsdf <- SpatialPointsDataFrame(spoints,df) # create a SpatialPointsDataFrame
plot(spointsdf,add=T,col=c('white','red'),pch=19,cex=0.5) # plot it on our map
text(121,24, 'TAIWAN', cex=1)

- Coordinates of these two university are not in spatial coordinate
format, which means it may run off while you change projection style of
your map
Fix position with coordinates
plot(wrld_simpl,xlim=c(-130,-60),ylim=c(45,80),col='#D2B48C',bg='lightblue')
coords <- matrix(c(-110,-102,-102,-110,-110,60,60,49,49,60),ncol=2)
l <- Line(coords)
ls <- Lines(list(l),ID="1")
sls <- SpatialLines(list(ls))
df <- data.frame(province="Saskatchewan") #create an empty dataframe to put lines into it.
sldf <- SpatialLinesDataFrame(sls,df)
plot(sldf,add=T,col='#3d2402', cex=2)
text(-114, 55, 'Saskatchewan', srt=90, cex=0.5)
text(-114, 63, 'CANADA', cex=1)

generate map from GADM data
library(raster)
TWN1 <- getData('GADM', country="TWN", level=0) # data Taiwan
## Warning in getData("GADM", country = "TWN", level = 0): getData will be removed in a future version of raster
## . Please use the geodata package instead
JPN <- getData('GADM', country="JPN", level=0) # data Japan
## Warning in getData("GADM", country = "JPN", level = 0): getData will be removed in a future version of raster
## . Please use the geodata package instead
class(TWN1) # those datasets are SpatialPolygonsDataFrame
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
par(mfrow = c(1, 2))
plot(TWN1,axes=T,bg=colors()[431],col='grey')
plot(JPN,axes=T,bg=colors()[431],col='grey')

Zoom on map
plot (TWN1, axes=T, xlim=c(121,122), ylim=c(24,25.5), bg=colors()[431],col='gray')

plot (TWN1, axes=T, xlim=c(121,122), ylim=c(24,25.5), bg=colors()[431],col='grey')

TWN2 <- getData('GADM', country="TWN", level=1)
## Warning in getData("GADM", country = "TWN", level = 1): getData will be removed in a future version of raster
## . Please use the geodata package instead
TWN2$NAME_1
## [1] "Fujian" "Kaohsiung" "New Taipei" "Taichung" "Tainan"
## [6] "Taipei" "Taiwan"
plot(TWN1,col="grey",xlim=c(119,122.5), ylim=c(21.5,25.5), bg=colors()[431], axes=T)
KAO <- TWN2[TWN2$NAME_1=="Kaohsiung",]
plot(KAO,col="grey 33",add=TRUE)

# base map
plot(TWN1,col="grey",xlim=c(121,122), ylim=c(24,25.5), bg=colors()[431], axes=T)
# adding spatial polygones
TAI <- TWN2[TWN2$NAME_1=="Taipei" | TWN2$NAME_1=="New Taipei" ,]
plot(TAI,col="black",add=TRUE)
# adding spatial points
coords <- matrix(cbind(lon=c(121.2,121.55,121.8),lat=c(25,25.19,24.5)),ncol=2)
coords <- coordinates(coords)
spoints <- SpatialPoints(coords)
df <- data.frame(location=c("City 1","City 2","City 3"),pop=c(138644,390095,34562))
spointsdf <- SpatialPointsDataFrame(spoints,df)
scalefactor <- sqrt(spointsdf$pop)/sqrt(max(spointsdf$pop))
plot(spointsdf,add=TRUE,col='white',pch=1,cex=scalefactor*3,lwd=2)
# adding a location of NTU (not spatial point here)
points(121.537290,25.021335, type="p", pch=18, col='white', cex=1.5)
# adding text
text(121.53,24.921335,"NTU", col='white', font=2)
# adding scale
maps::map.scale(x=120, y=25.4)

ggplot2
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6 ✔ purrr 0.3.4
## ✔ tibble 3.1.8 ✔ dplyr 1.0.10
## ✔ tidyr 1.2.1 ✔ stringr 1.4.1
## ✔ readr 2.1.3 ✔ forcats 0.5.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::extract() masks raster::extract()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ dplyr::select() masks raster::select()
library(rnaturalearth)
library(rnaturalearthdata)
library(sf)
## Linking to GEOS 3.9.1, GDAL 3.4.3, PROJ 7.2.1; sf_use_s2() is TRUE
theme_set(theme_bw())
world <- ne_countries(scale = "medium", returnclass = "sf")
class(world)
## [1] "sf" "data.frame"
ggplot(data = world) +
geom_sf() +
coord_sf(expand = FALSE) +
xlab("Longitude") + ylab("Latitude") +
ggtitle("World map", subtitle = paste0("(", length(unique(world$name)), " countries)")) +
geom_sf(color = "black", fill = "lightgreen")
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

dev.off()
## null device
## 1
GBIF data
library(rgbif)
library(mapr)
gbif.res <- occ_search(scientificName = "Urocissa caerulea", limit = 1200)
map_ggplot(gbif.res) +
coord_sf(xlim = c(120, 123), ylim = c(21, 26), expand = FALSE)

Bathymetric map with marmap
library(marmap)
## Registered S3 methods overwritten by 'adehabitatMA':
## method from
## print.SpatialPixelsDataFrame sp
## print.SpatialPixels sp
##
## 載入套件:'marmap'
## 下列物件被遮斷自 'package:raster':
##
## as.raster
## 下列物件被遮斷自 'package:grDevices':
##
## as.raster
# query
TW.bathy <- getNOAA.bathy(lon1=118,lon2=124, lat1=21,lat2=26,resolution=1) # don't put too wide / resolution: 1
## Querying NOAA database ...
## This may take seconds to minutes, depending on grid size
## Building bathy matrix ...
# define palette
blues <- colorRampPalette(c("darkblue", "cyan"))
greys <- colorRampPalette(c(grey(0.4),grey(0.99)))
# make the plot, step: difine line where to exist (e.g. draw lines per 2000 meter between -6000 to -1000 meters)
plot.bathy(TW.bathy,
image=T,
deepest.isobath=c(-6000,-120,-30,0),
shallowest.isobath=c(-1000,-60,0,0),
step=c(2000,60,30,0),
lwd=c(0.3,1,1,2),
lty=c(1,1,1,1),
col=c("grey","black","black","black"),
drawlabels=c(T,T,T,F),
bpal = list(c(0,max(TW.bathy),greys(100)),c(min(TW.bathy),0,blues(100))),
land=T, xaxs="i"
)

Leaflet
library(leaflet)
library(rgdal)
## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
##
## rgdal: version: 1.5-32, (SVN revision 1176)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.4.3, released 2022/04/22
## Path to GDAL shared files: C:/Users/ASUS/AppData/Local/R/win-library/4.2/sf/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/ASUS/AppData/Local/R/win-library/4.2/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.5-0
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
taiwan <- readOGR('D:/R_WD/Git linked/2022RisFUN-Mo/Mapping/mapdata202209220943/COUNTY_MOI_1090820.shp', use_iconv=TRUE, encoding='UTF-8')
## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS =
## dumpSRS, : Discarded datum Taiwan_Datum_1997 in Proj4 definition: +proj=longlat
## +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs
## OGR data source with driver: ESRI Shapefile
## Source: "D:\R_WD\Git linked\2022RisFUN-Mo\Mapping\mapdata202209220943\COUNTY_MOI_1090820.shp", layer: "COUNTY_MOI_1090820"
## with 22 features
## It has 4 fields
FRE <- paste(sep = "<br/>",
"<b><a href='https://www.dipintothereef.com/'>FRELAb TAIWAN</a></b>",
"Functional Reef Ecology Lab",
"Institute of Oceanography, NTU")
# Change style by addtiles
leaflet(taiwan) %>%
addPolygons(weight=0.5) %>%
addTiles(group="Kort") %>%
addPopups(121.53725, 25.021252, FRE, options = popupOptions(closeButton = FALSE))